//
// Alain Boucher - 19/03/2008
//
// Fonctions utilitaires pour la transformée de Fourier
// TP du cours de Traitement d'images
//
// Note : Les fonctions ci-dessous ne contiennent pas (volontairement)
//        de commentaires dans le code. Les seuls commentaires sont dans 
//        la facon d'utiliser les fonctions. Modifier les parametres et
//        observer les resultats pour bien comprendre le fonctionnement.
//

#include "ImageGray.h"
#include "fft2d_util.h"

//******   Structure d'un point *******//
// permet de manilpuler les points

using namespace Torch;


	struct Point
	{
	int abs;
        int ord;
        Point(int x=0, int y=0)
        {
     	 abs = x;
   	 ord = y;
         return ;
	} ;
        Point &operator=(const Point &P)
	{
    	 abs = P.abs;
   	 ord = P.ord;
   	 return *this;
	};
	};

//
// Tab1D_to_Image2D
//
// Cette fonction transforme un tableau 1D en une image 2D
//
// Entrees
// tab : tab 1D de taille width*height
// improve : rehausse le contenu de l'image par une transformation logarithmique
// swap : echange les quadrants pour mettre l'origine au centre de l'image
//
// Sortie
// une image ImageGray de taille width x height
//

ImageGray *Tab1D_to_Image2D (real *tab, int width, int height, bool improve, bool swap)
{
	int size = width * height;

	real *tab2 = new real [size];
	for(int i = 0 ; i < size ; i++)
		tab2[i] = fabs(tab[i]);
	
	if (improve)
		for(int i = 0 ; i < size ; i++)
			tab2[i] = log ( 1.0 + fabs (tab2[i]) );

	real min = 0.;
	real max = 0.;
	for(int i = 0 ; i < size ; i++)
	{
		if ((tab2[i] < min) || (i==0))  min = tab2[i];
		if ((tab2[i] > max) || (i==0))  max = tab2[i];
	}
	real maxmin = max - min;

	if (swap) SwapTabImage (tab2, width, height);

        ImageGray *image_tab = new ImageGray (width, height);

	for(int i = 0 ; i < size ; i++)
	{
		real val = ((tab2[i] - min) / maxmin) * 255.;
		if (val<0.)   val = 0.;
		if (val>255.)   val = 255.;
		image_tab->pixmap[i] = (unsigned char)val;
	}

	delete[] tab2;
	return image_tab;
}

//
// SwapTabImage
//
// Echange les quadrants d'une image pour mettre l'origine au centre (sinon en (0,0))
//
// 	AB  ->  DC
// 	CD      BA
//
// tab : tab 1D de taille width*height  (tableau modifie en sortie)
//
// Note : Pour les besoins du TP, vous n'avez pas besoin d'appeler cette fonction.
//        Cette fonction est automatiquement appelee dans Tab1D_to_Image2D.
//

void SwapTabImage (real *tab, int width, int height)
{
	int width2 = width/2;
	int height2 = height/2;

	real val = 0.;
	int ind1=0, ind2=0;
	for (int j=0; j<height2; j++)
	{
		for (int i=0; i<width2; i++)
		{
			ind1 = (j*width) + i;
			ind2 = ((j+height2) * width) + i + width2;
			val = tab[ind1];
			tab[ind1] = tab[ind2];
			tab[ind2] = val;
		}
		for (int i=width2; i<width; i++)
		{
			ind1 = (j*width) + i;
			ind2 = ((j+height2) * width) + i - width2;
			val = tab[ind1];
			tab[ind1] = tab[ind2];
			tab[ind2] = val;
		}
	}

	return;
}


//
// DefinirFiltre1
//
// Definit un filtre a appliquer sur une image frequentielle (apres Transformee de Fourier)
// Ce filtre est additif (addition pixel a pixel sur les images reelle et imaginaire)
//
// Entrees
// int width, height : taille de l'image a filtrer
// float facteur : parametre du filtre : 0 < facteur < 1
//
// Note : On définit ci-dessous un filtre binaire (0 ou 1). Pour accentuer les
// effets de ce filtre (et mieux comprendre), il est suggérer de remplacer les
// 1 par la valeur maximale de l'image lors du filtrage.
//
// Sortie
// un tableau 1D de taille width*height contenant le filtre
//

real *DefinirFiltre1 (int width, int height, float facteur)
{
	int size = width * height;

	real *filtre = new real [size];
	for (int i=0; i<size; i++) filtre[i] = 0;

	int facty = int( float(height/2) * facteur );
	int factx = int( float(width/2)  * facteur );

	filtre[(facty*width) + factx] = 1.;
	filtre[((height-facty)*width) + width - factx] = 1.;

	return filtre;
}

//
// AppliquerFiltre1
//
// Filtrage additif du tableau 1D fft par le filtre filtre1
//
// Note : Pour accentuer le filtrage, on remplace les valeurs 1 par le maximum de l'image.
//
// En sortie, le tableau fft est modifie

//

void AppliquerFiltre1 (real *fft, real *filtre1, int size)
{
	real max=0;

	for (int i=0; i<size; i++)
	   if ((fft[i]>max) || (i==0))  max = fft[i];

	for (int i=0; i<size; i++)
		if (filtre1[i] > 0.) fft[i] = max;

	return;
}

//
// AppliquerBackFiltre1
//
// Recheche des positions des points du tableau 1D fft qui on été préalablement modifiés par le filtre filtre1 
//
// Note : Pour revenir à l'image de départ, théoriquement il suffit de mettre tous ces points à 0.
//
// En sortie, un tableau contenant les positions de tous les points qui ont été modifiés par le filtres
//

int *AppliquerBackFiltre1 (real *fft, real *filtre1, int size)
{

        int *pos=new int [size];
       
        int j=0; 
	for (int i=0; i<size; i++)
	 if (filtre1[i] > 0.) 
          {
           pos[j] = i;
           j++;
          }	
	 pos[j] =-1;


	return pos;
}

//
// DefinirFiltre2
//
// Definit un filtre a appliquer sur une image frequentielle (apres Transformee de Fourier)
// Ce filtre est multiplicatif (multiplication pixel a pixel sur les images reelle et imaginaire)
//
// Entrees
// int width, height : taille de l'image a filtrer
// float facteur : parametre du filtre : 0 < facteur < 1
// bool inverse : parametre du filtre : true | false
//
// Note : On définit ci-dessous un filtre binaire (0 ou 1). 
//
// Sortie
// un tableau 1D de taille width*height contenant le filtre
//

real *DefinirFiltre2 (int width, int height, float facteur, bool inverse)
{
	real val1 = 0.;
	real val2 = 1.;
	if (inverse)
	{
		val1 = 1.;
		val2 = 0.;
	}

	int size = width * height;
	real *filtre = new real [size];
	for (int i=0; i<size; i++) filtre[i] = val2;

	int height2 = height/2;
	int width2  = width/2;

	int radius = int(float(height2)*facteur);
	if (width2 < height2)
		radius = int(float(width2)*facteur);

	int radius2 = radius*radius;

	for (int y=0; y<=radius; y++)
		for (int x=0; x<=radius; x++)
		{
			if (((y*y) + (x*x)) <= radius2)
			{
				filtre[(y*width) + x] = val1;
				if ((y>0) && (x>0)) filtre[((height-y)*width) + width-x] = val1;
				if (y>0)            filtre[((height-y)*width) + x] = val1;
				if (x>0)            filtre[(y*width) + width-x] = val1;
			}
		}

   return filtre;
}


//
// AppliquerFiltre2
//
// Filtrage multiplicatif du tableau 1D fft par le filtre filtre2
//
// En sortie, le tableau fft est modifie
//

void AppliquerFiltre2 (real *fft, real *filtre2, int size)
{
	for (int i=0; i<size; i++)
		if (filtre2[i] < 0.1) fft[i] = 0;

	return;
}


// noiseRemove
//
// Noircir tout les points du tableau fft qui ont été modifié par le filtre 1 : ces points sont déterminé par le tableau pos.
// Nous allons donc définir un contour sur chaque points et nous allons noicir toutes ces zone: exactement comme dans l'application du filtre 2
// Note: Etant donné que la convertion de l'image avec FFT inverse (cropGray) me générer pas deux points seulement je vais donc déternimer une ensemble
// de points sur la diagonale que je vais également noircir. 
//    
// En sortie, le tableau fftr ou ffti est modifie
//
real *noiseRemove (real *tab, int width, int height, bool improve, bool swap,int *pos)
{

  int size = width * height;
  real *filtre = new real [size];
  Point *Points;
  Points = new Point[8];
  int posy;	
  int posx;

 int k=0;
 int i=0;

// Calcul des coordonnées des deux points principales

  for(i = 0 ; k < 2 ; i++)
   {
     if(Appartient(pos,i))
     {	
      	posy = (int)(i/width);	
	posx = (int)(i%width);
	Point Ptaux(posx,posy);
        Points[k] = Ptaux;
        k++;
     }
   }

// Calcul des coordonnées des autres points générés par la procedure inverse

   int difx=  Points[0].abs;
   int dify=  Points[0].ord;
   int n=(int)(width-difx)/difx;
   n+=(n%2);  
   int max=40;  
   k=1;
   bool check=true;

  while(k<n && k<max) 
  {
     k++;
    if(k%2==0)
     {
	posy = Points[k-2].ord+dify;	
	posx = Points[k-2].abs+difx;
       
	if(true)
        {
	Point Ptaux1(posx,posy);
	Points[k] = Ptaux1;
        }
        else
        {
        check=false;
        break;
        }
     }
     else
     { 
     	posy = Points[k-2].ord-dify;	
	posx = Points[k-2].abs-difx;

	if(true)
        {
	Point Ptaux2(posx,posy);
	Points[k] = Ptaux2;  
        }
        else
        {
        check=false;
        break;
        }

     }  
   }

// Application du filtre sur ces points exacetement comme avec le filtre 2 mais avec pour facteur dépendant du nombre de points: pour qu'il n'y ai pas
// chevauchement entre zones.

   for(int i = 0 ; i < k ; i++)
   { 
        real fact=0.2/(n-1); 
        filtre=filtrePoint(width,height,fact,Points[i].abs,Points[i].ord); 
        AppliquerFiltre2 (tab,filtre,size);
 
   }
   return tab;
}

// Filtre les allentours d'un point
// Pour un point donné il filtre la table fft pour mettre tous les contours du point à 0.
// Il fonctionne comme le filtre 2, avec un facteur en entrée.

real *filtrePoint (int width, int height, float facteur,int posx,int posy)
{

	real val1 = 0.;
	real val2 = 1.;

	int size = width * height;
	real *filtre = new real [size];
	for (int i=0; i<size; i++) filtre[i] = val2;

	int height2 = height/2;
	int width2  = width/2;

	int radius = int(float(height2)*facteur);
	if (width2 < height2)
		radius = int(float(width2)*facteur);

	int radius2 = radius*radius;

	for (int y=0; y<=radius; y++)
		for (int x=0; x<=radius; x++)
		{
			if (((y*y) + (x*x)) <= radius2)
			{
		          filtre[((y+posy)*width) + x + posx] = val1;
			  filtre[((y+posy)*width) - x + posx] = val1;
			  filtre[((posy-y)*width) + x + posx] = val1;
			  filtre[((posy-y)*width) - x + posx] = val1;

			}
		} 
return filtre;
}

// verification d'appartenance

bool Appartient(int *pos,int val)
{
int i=0;
while(pos[i])
{
 if(pos[i]==val) return true;
 i++;    
}
return false;
}
// fin de fft2d_util.cc



